數值方法入門：以單擺為例理解 Euler、RK2（midpoint）、RK4（classical）

關鍵信息密度高、一步到位。你之後回看，應該可以即刻上手。

0. 問題背景與符號

考慮常微分方程（ODE）初值問題：

\dot{y} = f (t, y), y (t_{0}) = y_{0},

$f$ $t_n\to t_{n+1}=t_n+h$ $y_{n+1}$ 。

1. 從 Taylor 展開到 Euler／RK

以 Taylor 展開：

y (t + h) = y (t) + h f (t, y) + \frac{h^{2}}{2} f^{'} (t, y) + O (h^{3}) .

Euler（顯式）：只保留第一項
$y_{n + 1} = y_{n} + h f (t_{n}, y_{n})$
長方形積分 $O(h^2)$ $O(h)$ 。
RK2（midpoint）：取中點斜率
$\begin{aligned} k_{1} & = f (t_{n}, y_{n}), \\ k_{2} & = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}), \\ y_{n + 1} = y_{n} + h k_{2} \end{aligned}$
梯形／中點積分 $O(h^2)$ 。
RK4（classical）：四點加權平均
$\begin{aligned} k_{1} & = f (t_{n}, y_{n}), \\ k_{2} & = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}), \\ k_{3} & = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}), \\ k_{4} & = f (t_{n} + h, y_{n} + h k_{3}), \\ y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) \end{aligned}$
多點平均斜率 $O(h^4)$ 。

2. 單擺（Pendulum）作為非線性測試牀

模型（無阻尼）

\begin{matrix} \ddot{θ} + \frac{g}{L} \sin θ = 0 \Rightarrow y = [\begin{matrix} θ \\ ω \end{matrix}], f (y) = [\begin{matrix} ω \\ - \frac{g}{L} \sin θ \end{matrix}] . \end{matrix}

2.1 矩陣（向量）更新式（Euler／RK2／RK4）

Euler：

\begin{matrix} y_{n + 1}^{(Euler)} = [\begin{matrix} θ_{n} + h ω_{n} \\ ω_{n} - \frac{g h}{L} \sin θ_{n} \end{matrix}] \end{matrix}

RK2（Midpoint）：

\begin{matrix} y_{n + 1}^{(RK 2)} = [\begin{matrix} θ_{n} + h ω_{n} - \frac{g h^{2}}{2 L} \sin θ_{n} \\ ω_{n} - \frac{g h}{L} \sin (θ_{n} + \frac{h ω_{n}}{2}) \end{matrix}] \end{matrix}

RK4（Classical） $k_1,\dots,k_4$ 的矩陣式寫法：

\begin{aligned} k_{1} & = [\begin{array}{c} ω_{n} \\ - \frac{g}{L} \sin θ_{n} \end{array}], k_{2} = [\begin{array}{c} ω_{n} + \frac{h}{2} k_{1, ω} \\ - \frac{g}{L} \sin (θ_{n} + \frac{h}{2} k_{1, θ}) \end{array}], \\ k_{3} & = [\begin{array}{c} ω_{n} + \frac{h}{2} k_{2, ω} \\ - \frac{g}{L} \sin (θ_{n} + \frac{h}{2} k_{2, θ}) \end{array}], k_{4} = [\begin{array}{c} ω_{n} + h k_{3, ω} \\ - \frac{g}{L} \sin (θ_{n} + h k_{3, θ}) \end{array}], \end{aligned}

2.2 真解（作對照）

Jacobi 橢圓函數 $\sin\theta\!\approx\!\theta$ $\theta(t)\approx\theta_0\cos(\sqrt{g/L}\,t)$ 。

由上圖可見，如果h的值不夠細，即計算精度不夠高，在Euler算法上會嚴重偏移，但在RK4上還保持良好。

3. Gravity ODE：何時需要 RK？何時不需要？

恆定重力： $\dot x=v,\ \dot v=a=(0,0,-g)$ 解析解：

x (t) = x_{0} + v_{0} t + \frac{1}{2} a t^{2}, v (t) = v_{0} + a t .

→ 這時 $O(h)$ 誤差；RK2 已精確；RK4/Adaptive RK 沒必要。 需要 RK 的情形：

$v$ ）：非線性、無解析解
$g(r)=GM/r^2$ ）：非線性
N-body 多體引力：無解析解，只能數值法（甚至要保結構積分器）

4. 輕量可重用工具：固定步長 Euler／RK2／RK4（Python）

我幫你寫咗通用整合器模組（最小但實用）：

下載模組： ode_integrators.py
介面極簡：


xxxxxxxxxx
2
1
from ode_integrators import integrate
2
ts, ys = integrate(f, t0=0.0, y0=np.array([...]), t1=10.0, h=0.01, method="rk4")

method 可選 "euler" | "rk2" | "rk4"；固定步長 h；回傳等距 ts、配對 ys。

4.1 範例 A：恆定重力（6 維）


xxxxxxxxxx
12
1
import numpy as np
2
from ode_integrators import integrate
3

4
g = 9.81
5
def grav_rhs(t, y):
6
    # y = [x, y, z, vx, vy, vz]
7
    x, y_, z, vx, vy, vz = y
8
    ax, ay, az = 0.0, 0.0, -g
9
    return np.array([vx, vy, vz, ax, ay, az], dtype=float)
10

11
y0 = np.array([0, 0, 0,   10, 5, 20], dtype=float)
12
ts, ys = integrate(grav_rhs, t0=0.0, y0=y0, t1=2.0, h=0.1, method="rk2")

4.2 範例 B：單擺（2 維）


xxxxxxxxxx
13
1
import numpy as np
2
from ode_integrators import integrate
3

4
g, L = 9.81, 1.0
5
def pend_rhs(t, y):
6
    th, om = y
7
    return np.array([om, -(g/L)*np.sin(th)], dtype=float)
8

9
y0 = np.array([np.deg2rad(45.0), 0.0])
10
ts, ys = integrate(pend_rhs, 0.0, y0, 0.5, 0.05, method="rk4")
11

12
theta_deg = np.rad2deg(ys[:,0])
13
omega_deg_s = np.rad2deg(ys[:,1])

想要事件偵測（撞地、過零交叉）、能量監控、或自動步長？請看下一節 SciPy。

5. 現成高階解算器：SciPy `solve_ivp`（自動步長／誤差控制）

如可用 SciPy，建議一般專案直接用 solve_ivp：


xxxxxxxxxx
24
1
import numpy as np
2
from scipy.integrate import solve_ivp
3

4
# 單擺例子
5
g, L = 9.81, 1.0
6
def f(t, y):
7
    # y = [theta, omega]
8
    th, om = y
9
    return np.array([om, -(g/L)*np.sin(th)], dtype=float)
10

11
t_span = (0.0, 10.0)
12
y0 = np.array([np.deg2rad(45.0), 0.0], dtype=float)
13

14
# 常用顯式自適應：RK45（Dormand–Prince）
15
sol = solve_ivp(f, t_span, y0, method="RK45",
16
                rtol=1e-6, atol=1e-9, dense_output=True)
17

18
# 固定時間網格取樣
19
t_eval = np.linspace(0, 10, 501)
20
sol = solve_ivp(f, t_span, y0, method="RK45",
21
                t_eval=t_eval, rtol=1e-6, atol=1e-9)
22

23
ts = sol.t           # (N,)
24
ys = sol.y.T         # (N, d)  注意：SciPy 的 y 是 (d, N)

方法選擇指引：

RK45：大多數非剛性問題的默認選擇
DOP853：更高階顯式 RK（步數更省）
Radau / BDF：剛性問題（隱式，穩定域大）
rtol/atol：誤差控制旋鈕；越嚴格越慢但更準

Gravity 的建議

恆定重力：closed-form 或我上面的 rk2 已經夠；
變重力場／空阻／N-body：建議 RK45/DOP853；剛性時 Radau/BDF。

6. 常見坑位 Checklist

步長是時間（秒），唔係角度或空間單位
單位一致（deg vs rad、m vs cm、N vs kN）
能量漂移：Euler 會嚴重漂移；RK2 好很多；RK4 更穩
剛性系統：顯式發散 → 試 Radau/BDF
事件偵測：solve_ivp 的 events 可做（撞牆、過零）
多體問題：注意數值穩定性與守恆（考慮 symplectic integrator）

附：可直接貼入你筆記的 SymPy 檢核程式（等價驗證＆輸出 LaTeX）


xxxxxxxxxx
52
1
# SymPy: Pendulum ODE updates in matrix form + RK4 equivalence check
2

3
import sympy as sp
4

5
theta_n, omega_n, h, g, L = sp.symbols('theta_n omega_n h g L', real=True)
6
y_n = sp.Matrix([theta_n, omega_n])
7

8
def f(y):
9
    th, om = y
10
    return sp.Matrix([om, -(g/L) * sp.sin(th)])
11

12
# Euler
13
y_euler = sp.simplify(y_n + h * f(y_n))
14

15
# RK2 (midpoint)
16
k1_rk2 = f(y_n)
17
k2_rk2 = f(y_n + (h/2) * k1_rk2)
18
y_rk2 = sp.simplify(y_n + h * k2_rk2)
19

20
# RK4 (teaching form)
21
k1 = f(y_n)
22
k2 = f(y_n + (h/2) * k1)
23
k3 = f(y_n + (h/2) * k2)
24
k4 = f(y_n + h * k3)
25
y_rk4_A = sp.simplify(y_n + (h/6) * (k1 + 2*k2 + 2*k3 + k4))
26

27
# RK4 (long-form stages)
28
theta_1 = theta_n + (h/2)*k1[0];   omega_1 = omega_n + (h/2)*k1[1]
29
k2_alt  = sp.Matrix([omega_1, -(g/L)*sp.sin(theta_1)])
30

31
theta_2 = theta_n + (h/2)*k2_alt[0]; omega_2 = omega_n + (h/2)*k2_alt[1]
32
k3_alt  = sp.Matrix([omega_2, -(g/L)*sp.sin(theta_2)])
33

34
theta_3 = theta_n + h*k3_alt[0];    omega_3 = omega_n + h*k3_alt[1]
35
k4_alt  = sp.Matrix([omega_3, -(g/L)*sp.sin(theta_3)])
36

37
y_rk4_B = sp.simplify(y_n + (h/6) * (k1 + 2*k2_alt + 2*k3_alt + k4_alt))
38

39
# Symbolic & numeric checks
40
print("Symbolic diff:")
41
sp.pprint(sp.simplify(y_rk4_A - y_rk4_B))
42

43
vals = {theta_n: 0.7, omega_n: -0.3, h: 0.05, g: 9.81, L: 1.2}
44
print("\nNumeric spot-check:", vals)
45
print("y_rk4_A =", sp.N(y_rk4_A.subs(vals)))
46
print("y_rk4_B =", sp.N(y_rk4_B.subs(vals)))
47
print("difference =", sp.N((y_rk4_A - y_rk4_B).subs(vals)))
48

49
# LaTeX blocks
50
print("\nLaTeX Euler:\n", sp.latex(y_euler))
51
print("\nLaTeX RK2:\n",   sp.latex(y_rk2))
52
print("\nLaTeX RK4:\n",   sp.latex(y_rk4_A))